function sa2(N)

  close all
  

  function resout = outpt(t, y, flag)
    disp(t)
    resout = 0;
  end
  
  %time_steps = [0:N-1]./(N-1)*1.0;
  z = chebpoints(N, 0, 0.999);

  function resSa2a = Sa2a(z)
    if (z < 0.001)
      resSa2a = (4/3)+(-3/2).*z+(5/3).*z.^3+(-12/7).*z.^4+(16/9).*z.^6+(-9/5).*z.^7+(11/6).*z.^9+(-24/13).*z.^10+(28/15).*z.^12+(-15/8).*z.^13+(17/9).*z.^15+(-36/19).*z.^16+(40/21).*z.^18+(-21/11).*z.^19; % /z^2
    else
      resSa2a = (1/3).*z.^(-3).*(3.^(1/2).*pi+(-6).*z.*(1+z).*(1+z+z.^2).^(-1)+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)); % / z^2
    end
  end
  function resSa2b = Sa2b(z)
    if (z < 0.001)
      resSa2b = (-1/3)+(3/8).*z+(-5/12).*z.^3+(3/7).*z.^4+(-4/9).*z.^6+(9/20).*z.^7+(-11/24).*z.^9+(6/13).*z.^10+(-7/15).*z.^12+(15/32).*z.^13+(-17/36).*z.^15+(9/19).*z.^16+(-10/21).*z.^18+(21/44).*z.^19;
    else
      resSa2b = z.^(-2).*((1/2)+(-1/12).*z.^(-1).*(3.^(1/2).*pi+6.*z.^3.*(1+z+z.^2).^(-1)+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)));
    end
  end
  function resSa2c = Sa2c(z)
    if (z < 0.001)
      resSa2c = 1+(-9/8).*z+(5/4).*z.^3+(-9/7).*z.^4+(4/3).*z.^6+(-27/20).*z.^7+(11/8).*z.^9+(-18/13).*z.^10+(7/5).*z.^12+(-45/32).*z.^13+(17/12).*z.^15+(-27/19).*z.^16+(10/7).*z.^18+(-63/44).*z.^19;
    else
      resSa2c = z.^(-2).*((-1/2)+(1/4).*(2.*((-2)+(-2).*z+z.^2).*(1+z+z.^2).^(-1)+z.^(-1).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2))));
    end
  end

  function resM = M(z)
  % pulled out z^3
    resM = zeros(2,2);
    resM(1,1) = 1;
    resM(2,1) = (-4+7.0.*z^3);
    resM(2,2) = z.*(z.^3-1);
  end

  function res = deq(z, u, t)
    a   = u(1);
    az  = u(2);
    res = zeros(2,1);
    res(1) = az;
    if (t == 1)
      res(2) = 2*Sa2a(z) - 9.0.*z.^2.*a;
    elseif (t == 2)
      res(2) = 2*Sa2b(z) - 9.0.*z.^2.*a + 2.0.*z;
    else 
      res(2) = 2*Sa2c(z) - 9.0.*z.^2.*a - 2.0.*z;
    end
  end

  reltol = 1e-8;


  % a
  uin = [(2+pi./3./sqrt(3)-log(3))./3; -4./6];
  opts = odeset('RelTol', reltol, 'AbsTol', 1e-14, 'MStateDependence', 'none', 'Mass', @M, 'OutputFcn', @outpt);
  [za, uout] = ode15s(@(z, u) deq(z, u, 1), z, uin, opts);
  a = uout(:, 1);
  na2a = a.*za.^3;
  na2az = uout(:, 2);

  % b
  uin = [(1/3)+(1/108).*((-1).*3.^(1/2).*pi+9.*((-2)+log(3))); -1./4];
  opts = odeset('RelTol', reltol, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @M);
  [zb, uout] = ode15s(@(z, u) deq(z, u, 2), z, uin, opts);
  a = uout(:, 1);
  na2b = a.*zb.^3;
  na2bz = uout(:, 2);

  % c
  uin = [(-1/3)+(1/36).*(18+3.^(1/2).*pi+(-9).*log(3));1./4];
  opts = odeset('RelTol', reltol, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @M);
  [zc, uout] = ode15s(@(z, u) deq(z, u, 3), z, uin, opts);
  a = uout(:, 1);
  na2c = a.*zc.^0;
  na2cz = uout(:, 2);

  na2c = SpectralGrad(na2c,'cheb2'); 

  figure
  plot(za, na2a, zb, na2b, zc, na2c);

  figure
  plot(za, na2az, zb, na2bz, zc, na2cz);
%plot(z, res);
  %uout(:,1)./log(1-tout)


end

